function [ys,check] = CGH_steadystate(ys,exe)
global M_ normalizationI fixed_costI 

% Read out parameters
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

%Initialize indicator
check = 0;
subsidy = 1;

load indicators

%%Initialize Using Analytical Steady State

siggma                  = 80*ezI + (1-ezI)*2;
theta_v                 = (1-siggma)/(1-1/psii);
V                       = 1;
E_t_V_tp1_1_minus_sigma = 1;
PHI                     = 1;
Pi                      = Pi_bar;
R                       = Pi_bar/betta;
Z_bar                   = 1;
Z                       = Z_bar;
a                       = a_bar;
sigma_a                 = sigma_a_bar;

u                       = 1;
u_n                     = 1;
Y                       = 1;
N_n                     = 0.19*homeprodI + (1-homeprodI)*1e-10; 
N_m                     = 0.33; 
K_m                     = 5.15*Y*capitalI + 1 - capitalI;
K_n                     = 6.75*Y*home_capitalI*capitalI + 1 - home_capitalI;

q_n                     = 1;
q                       = 1;
M                       = betta;
R_R                     = 1/betta;
I_m                     = delta_0*K_m;
I_n                     = delta_n_0*K_n*home_capitalI;
mu                      = theta_mu/(theta_mu-1);
Xi                      = 1/mu;
delta_1                 = 1/betta-1+delta_0;
delta_n_1               = 1/betta-1+delta_n_0;
R_K                     = delta_1;
Psi                     = (mu - 1)*fixed_costI;
alppha                  = R_K*K_m*(1/Xi)*((1+Psi)^-1);
PF_normalization        = normalizationI*(1/(Xi*N_m^(1-alppha)*K_m^alppha)) + (1-normalizationI);
W                       = ((1-alppha)*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/N_m;
D                       = Y - N_m*W - delta_0*K_m;
B                       = nu*K_m;
D_E                     = D-(1-1/R_R)*B;
P_E                     = betta/(1-betta)*D_E;
C_m                     = - I_n + - (P_E+1/R_R*nu*K_m)*capitalI + W*N_m + D_E + P_E*capitalI + nu*K_m;

if homeprodI==1
   alpha2P              = (K_n/W)*(delta_n_1/N_n)*(1+(K_n/W)*(delta_n_1/N_n))^-1;
else
   alpha2P              = 0;
end

C_n                     = ((u_n*K_n)^alpha2P)*(N_n^(1-alpha2P));
if homeprodI==1
    alpha1P                 = ((alpha2P/delta_n_1)*(C_n/K_n)*((C_n/C_m)^(b1P-1)))*((1 + (alpha2P/delta_n_1)*(C_n/K_n)*((C_n/C_m)^(b1P-1)))^-1);
else
    alpha1P                 = 1;
end
C                       = (alpha1P*(C_m^b1P) + (1-alpha1P)*(C_n^b1P))^(1/b1P);
L                       = 1 - N_m - N_n;
eta                     = ((W*alpha1P)*(L/C)*(C_m/C)^(b1P-1)+1)^-1;


options = optimset('Display','off','MaxIter', 5000,'MaxFunEvals',5000,'TolFun', 1.0000e-14,'TolX', 1.0000e-14);
x0 = [Y,alppha,eta,alpha1P,alpha2P];

[steadySol,fval]    = fsolve(@(x) CGH_steadystate_solve(x),x0,options);

Y       = steadySol(1);
alppha  = steadySol(2);
eta     = steadySol(3);
alpha1P = steadySol(4); 
alpha2P = steadySol(5); 


if max(abs(fval))>1e-8
    display('Inaccurate Steady State')
    keyboard
end

if min(steadySol)<-0.0001
    display('Negative Calibrated Value')
    keyboard
end


u                = 1;
u_n              = 1;
N_n              = 0.19*homeprodI + (1-homeprodI)*1e-10; 
N_m              = 0.33; 
K_m              = 5.15*Y*capitalI + 1 - capitalI;
K_n              = 6.75*Y*home_capitalI*capitalI + 1 - home_capitalI;

M                = betta;
mu               = theta_mu/(theta_mu-1);
Xi               = 1/(subsidy*mu);
Z                = 1;
u                = 1;
u_n              = 1;
q_n              = 1;
q                = 1;
delta_1          = 1/betta-1+delta_0;
delta_n_1        = 1/betta-1+delta_n_0;
R_R              = 1/betta;

if capitalI==1
    PF_normalization = normalizationI*(1/(Xi*N_m^(1-alppha)*K_m^alppha)) + (1-normalizationI);
    Psi              = fixed_costI*(1-Xi)*PF_normalization*N_m^(1-alppha)*K_m^alppha;
    I_m              = delta_0*K_m;
    I_n              = delta_n_0*K_n*home_capitalI;
    W                = ((1-alppha)*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/N_m;
    R_K              = (alppha*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/(u*K_m);
    D                = Y - N_m*W - delta_0*K_m;   
    B                = nu*K_m;
    D_E              = D-(1-1/R_R)*B;
    P_E              = betta/(1-betta)*D_E;
    C_m              = - I_n + - (P_E+1/R_R*nu*K_m) + W*N_m + D_E + P_E + nu*K_m;
else
    PF_normalization = normalizationI*(1/(Xi*N_m)) + (1-normalizationI);
    Psi              = fixed_costI*(1-Xi)*PF_normalization*N_m;
    I_m              = 0;
    I_n              = 0;
    W                = Xi*PF_normalization*(Z*N_m)/N_m;
    D                = Y - N_m*W;   
    R_K              = 1;
    B                = 0;
    D_E              = D;
    P_E              = 0;
    C_m              = - I_n + W*N_m;
end

if homeprodI==1
    if home_capitalI==1
       C_n              = ((u_n*K_n)^alpha2P)*(N_n^(1-alpha2P));
    else
       C_n              = N_n;    
    end
    C                = (alpha1P*(C_m^b1P) + (1-alpha1P)*(C_n^b1P))^(1/b1P) ;
else
    C_n = 1e-10;
    C                = C_m ;
end

L                = 1 - N_m - N_n;
utilC            = eta*(C^(eta-1))*(L^(1-eta));
utilL            = (1-eta)*(C^eta)*(L^-eta);
R_K              = (alppha*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/(u*K_m);
util             = C^eta*(L^(1-eta));
V_normalization  = ezI*((1-betta)/(util^((1-siggma)/theta_v))) + (1-ezI)*(1-betta)/(util^(1-siggma)/(1-siggma));


if abs(- C_m  + Y - I_m - I_n) >1e-7
    display('Resource Constraint Does Not Hold')
    keyboard
end

                      
sigma_z              = sigma_z_bar;  
PF_normalization     = log(PF_normalization);
PF_normalization_bar = PF_normalization;
a                    = 1;
Z                    = 1;
y_ss                 = Y; 

C_m_bar               = C_m;        
C_bar                 = C;        
C_n_bar               = C_n;                                                                                                                        
N_n_bar               = N_n;                                                                 
I_n_bar               = I_n;                                                                 
K_n_bar               = K_n;                                                                 
q_n_bar               = q_n;                                                                 
u_n_bar               = u_n;  
W_bar                 = W;
Lambda                = ezI*((V/util)^(1-(1-siggma)/theta_v)*utilC*alpha1P*(C_m/C)^(b1P-1)*a) +...
                        (1-ezI)*(util^-siggma)*utilC;
                    
alpha1_ex             = 0;

aux_sigma             = sigma_a_bar;


rho_sigma_alpha1     = rho_sigma_a;
sigma_alpha1         = sigma_a;
sigma_alpha1_bar     = sigma_a;
sigma_sigma_alpha1   = sigma_sigma_a;

load fig_number
if figure3==1 || figure4==1
    rho_sigma_a       = 0.74227; 
end
if figure5==1 
    rho_sigma_a    = 0.45*0.74227;
end

R_E=(D_E + P_E)/P_E;
E_R_E=R_E(+1);
E_R_E_squared=R_E^2;
cond_var_R_E=E_R_E_squared-E_R_E^2;
E_M_tp1=M;
E_R_E_risk_neutral=R_E/E_M_tp1;
E_M_tp1_R_E=M(+1)*R_E;
E_M_tp1_R_E_squared=M*R_E^2;
cond_var_R_E_risk_neutral=E_M_tp1_R_E_squared/E_M_tp1-(E_M_tp1_R_E/E_M_tp1)^2;

dNm = 0;
dNn = 0;
dL  = 0;

save steady y_ss


%% end model equations
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end

end
